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Renormalization Group Equations in integro-differential form describing the evo- 
lution of cascades or resumming logarithmic scaling violations have been known 
in quantum field theory for a long time. These equations have been traditionally 
solved by turning to Mellin moments, since in this space they become algebraic, 
x-space solutions are less known, but special asymptotic expansions exist which 
allow a fast numerical implementation of these equations. We illustrate how the 
equations can be solved using recursion relations in the ncxt-to-lcading order ap- 
proximation. 



1. Introduction 

^ Evolution equations describing the high energy behaviour of scattering 
amphtudes carry significant information on the factorization/ renormaUza- 
tion scale dependence of such amplitudes, and allow to link the behaviour 
of processes at a given energy scale to collisions taking place at another 
(usually much higher) scale. 

In QCD, the accepted theory of the strong interactions, confinement for- 
bids the detection of the fundamental states of the theory, such as quarks 
and gluons. However, asymptotic freedom allows to separate the pertur- 
bative dynamics at short distances from the non-perturbative one, due to 
confinement, through factorization theorems. The functions evolved by 
these equations are called parton distributions. We will very briefly intro- 
duce them below, and we will exploit the density matrix formulation of 
quantum mechanics as an analogy to illustrate the matter. 
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2. Nonlocal correlators, Wigner functions and all that 

The mathematical construct which is the closest to a parton distribution 
function (p.d.f ) q{x, Q) is a Wigner function. The analogy is, of course, 
limited. 

We recall that Wigncr's description of quantum mechanics via quasi- 
probabilities of phase space (x, p) 

/(x,p) = ^/ dy^*(x-|yMx+|y) (1) 

is fully equivalent to Schrodinger's formulation ^. 

Differently from Wigner functions, in a p.d.f a variable .x, called 
"Bjorkcn x" , takes now the role of the momentum "p" . Also, parton distri- 
butions are correlation functions of a special type, being defined just on the 
light-cone. In this sense they are not generic nonlocal correlators. This lim- 
itation, duo to the special nature of high c;nergy collisions in asymptotically 
free theories, sets the boundary of validity of the parton model approach 
to QCD. 

In QCD one starts by introducing, via arguments based on unitarity, 
the hadronic tensor, which is the key construct describing the collision 

^""^ / 0^A''"''^<PASA;PBSB\[J>'{(i),r{x)]\PASA;PBSB), (2) 

with Pa and Pg being the momenta of the colliding hadrons and Sa 
and Sb their covariant spins. The J's are electromagnetic currents. 

The distribution functions that emerge -at leading order- from this fac- 
torized picture arc correlation functions of non-local operators in configura- 
tion space. They are the quark-quark and the quark-antiquark correlators. 

Their expression simplifies in the axial gauge, in which the gauge link is 
removed by the gauge condition. For instance, the quark-quark correlator 
takes the form 

{^a/A)„0{P,S,k) = J ^e^'=-(P,5|4"^(0)V(»H^)|P5). (3) 

In (3) we have included the quark flavour index a and an index A for the 
hadron, as usual. Fields are not time ordered since they can be described 
by the good light cone components of ijj and by At, a transverse component 
of the gauge field A^, as discussed in ^. The non-perturbative information 
in a collision is carried by matrix elements of this type. 

Further considerations allow to show that the leading contributions to 
(3) come from the light-cone region. The leading expansion of the quark- 
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quark correlator then is of the form 




where we have used all the four vectors at our disposal (spin S, momentum 
P of the hadron) and introduced invariant amplitudes (parton distributions) 
/i, gi, hi, now expressed in terms of a scaling variable x (Bjorken x). 
is a light-cone four- vector (n^ = 0), approximately perpendicular to the 
hadron momentum. 

The definition of p.d.f.'s in (4) involves also an underlying physical scale 
Q {Q » Aqcd, with Aqcd being the scale of confinement), not appar- 
ent from that equation and characterizing the energy scale at which these 
matrix elements, summarized by (4), are defined. Truly: /i = fi{x,Q'^), 
hi — hi{x,Q'^) and so on. 

The role of the various renormalization group equations is to describe 
the perturbative change in these functions as the scale Q is raised (lowered). 
Each equation involves kernels {P{x)) of various types, of well known form, 
and asymptotic expansions of the solutions exist (see for instance ^ and the 
implementation given in ^). Here, however, we will illustrate an alternative 
method to solve these equations which is computationally very efficient. 

3. The equations 

Parton distributions, though fully identified by their operatorial form, are 
not currently calculable from first principles. In fact, the theory that they 
describe is a (nonlinear) gauge theory characterized by a large QCD cou- 
pling constant a at the scale at which they are usually introduced. On the 
other hand, the equations describing their evolution under the renormal- 
ization group are identified using perturbation theory 



Pi-){x,a{Q^))^Xi{x,Q^) 



(5) 



having defined 





with 
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(7) 
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for the non-singlet distributions and 

^ dQ^ \ G{x,Q^) ) 2-K \Pggix,Q^) Pgg{x,Q^)J \ G {x , Q"") 

(8) 

in the singlet sector. Here, G{x,Q'^) is the gluon density, while q{x,Q'^) is 
the quark density, a is the QCD coupling constant. 

Similar RG equations can be derived for the photon structure function 
q^{x,Q^), now with an inhomogeneus term included (aem is the QED cou- 
pling constant) 



Let's start from the latter equation. The ansatz for the solution of (9) is 
chosen of the form 

Q^) = ^(^±^ln-(^)+±^ln-(^)] (10) 

and the recursion relations for the fimctions An{x) and Bn{x) appearing 
in the expansion are obtained by comparing terms of the same order in 
a'' log^{a/ao), with k=0,l and n = 0, 1, .... We use a running QCD cou- 
pling at the desired perturbative order 

/^(") = -7r" -T^" (11) 



dlog(g2) ' 4n 167r2 

and ao = a{Qo), with Qo being the initial scale at which the evolution 
starts. In (11) (3o and /3i are the first two coefiicients of the QCD beta 
function. 

The recursion relations arc in leading order given by 

A^{x) = Ao{x) - K^°\x) - ^P^°\x) ^ Ao{x) 

Am+i{x) = Am{x) - ^P^'^Kx) (g) Am{x) m = 1, 2, 3... (12) 
Po 

and 

^i(^) = ^ {Am+i{x) - A^{x)) - 2^^-^ - ^P^'\x) ® Bm{x) - ^P^'\x) ® ^(o)(a;) 



'0 
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Po Po Po 

(13) 

at order a^. The intial condition is easily shown to be of the form 

with go (a;) defining the initial functional choice for the parton distribution 
at the lowest scale. It is possible to show ^ that the freedom in choosing the 
original values for Ao{x) and Bq{x) is not relevant at the order at which 
we are working (a^), as far as (14) is satisfied. 

In the case of other equations, such as eqs. (5), we get the recursion 
relations 

A„+i(:c) = -|-PW®A„(x) (15) 
Po 

and 

5„+i(a;) = -Bn{x) - (^An+iix)) - jV^''' ® M^) - ^P^°'> ® Bn{x) 
V4po / 47r/?o Po 



-Sn(a;)+(^P(°)®A„(X)) 

. 1 P(i) ^ A„(x) - -^p(o) ^ Bn{x), 
47r/?o Po 



(16) 



which are solved with the initial condition Bo{x) = 0. The initial conditions 
for the coefficients Ao{x) and Bo(x) are specified in a slightly different way 
from the photon case, with g(a;,(5o) now identified as the leading order 
ansatz for the initial distribution 

Ao{x)=S{l-x)(E)qix,Ql)=qo{x) (17) 

which also requires Bq{x) = 0, since we have to satisfy the boundary con- 
dition 

Ao{x) + QoSo(x) = qo{x). (18) 

Again, other boundary choices are possible for Aq{x) and Bo{x) as far 
as (18) is fuimUed. 
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4. Distributional Singularities and Finite Elements 

Once the recursion are given, it remains to be seen how to actually im- 
plement the method. In practice it is not so easy, but it is easier than in 
other methods ^. The codes can range from 1,000 lines for QCD to several 
thousands lines for supersymmctry. The advantage of these implementa- 
tions is that the codes can run in few minutes (2-3), compared to the much 
slower codes obtained before. The reason of such improvement is related 
to the use of analytical expressions (in a finite clement discretization of the 
integrals) which drastically reduce the computational time required to ac- 
tually perform the recursive integrations. There are some important points, 
however, to keep into account. The kernels P{x) are defined, in fact, in a 
distributional sense and arc plagued with artificial numerical singularities 
the most critical ones being tied to plus ( "+" ) distributions and defined as 

Jo Jo 1-2/ 

A simple trick to eliminate this singular behaviour is to use the identity 

' .«/(.). r^^^M-^+/(.)iog(i-.) (20) 



(1 -x)+ J^ y y- 

and proceed with a finite clement discretization of the resulting integral. 

We briefly recall the numerical strategy employed in this analysis. We 
define P{x) = xP{x) and A{x) = xA{x). We also define the convolution 
product 

J(x) ^ [ f (^) P (^) Aiy). (21) 

The integration interval in y at any fixed x-value is partitioned in an array 
of increasing points ordered from left to right {xq, Xi,X2, Xn, Xn+i) with 
xq = X and Xn+i = 1 being the upper edge of the integration region. One 
constructs a rescaled array {x, x/x^, ...,x/x2tx/xi, 1). We define Sj = x/xi, 
and Sn+i = X < Sn < Sn-i < ■■■si < So = 1. Wc get 

^«=Erf(;)K;)^*') 

At this point we introduce the linear interpolation 

A{y) =(l- ) A{x,) + A{x,+,) (23) 

\ Xi-\-\ Xi J Xi-\-\ X% 

and perform the integration on each subinterval with a change of variable 
y— > x/y and replace the integral J{x) with its discrete approximation 
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Jn{x) to get 

J^{x)=A{xo)-^ f^P{y){y-sr) 

^ Si- Si+1 Js,+, y 

j_i *» Jsi y 

(24) 

Introducing the coefficients W{x,x) and W{xi,x), the integral is cast in 
the form 

n 

Jn{x) = W{x, x)A{x) + W{x^,x)A{xi) (25) 

where 

W{x,x) = -^ f'^{y^si)P{y), 

1 - Si J SI V 

W{Xi,x)= r ^{y-Si+i)P{y) 

Si - Si+1 Js,+i y 

'-^ r^^{y-s,.i)P{y). 

Si-i - Si J^. y 

(26) 

The results of the integration of the recursion relations can be given 
in analytical form, with obvious care in handling the "+" distributions. 
We have solved by this method all the leading twist evolution equations of 
QCD to higher order. Application of the method to supersymmetry has 
also been illustrated ^. 

5. Applications 

As a simple illustration of the method we have included two figures 
which illustrate the variation in shape of the evolved functions with 
a varying final evolution scale. Scaling violations are usually quite 
small, however they are very important both to uncover important new 
physics and for precision studies. Good algorithms are always welcomed. 
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Figure 1. Plot of the supersymmetric evolution of the up and down quarks for a final 
scale of 100 GeV in a specific model. 
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Figure 2. Density of gluons, obtained by the recursion method for the QCD (reg.) and 
Supcrsynimetric QCD (SQCD) evolution. To vary here is the mass of the supersymmetric 
partner, the gluino, defining a matching scale in the QCD/SQCD evolution. 



